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ABSTRACT 

We numerically model the interaction between an expanding fireball and a stationary 
external medium whose density is either homogeneous or varies with distance as a power-law. 
The evolution is followed until most of the fireball kinetic energy is converted into internal 
energy. The density, pressure and flow Lorentz factor profiles are shown at different stages, 
including shock and rarefaction wave reflections, for a fireball of initial bulk Lorentz factor 
r = 100, both in the adiabatic and non-adiabatic (radiative) regimes. For cooling times shorter 
than the dynamic time, bolometric light-curves are computed for values of T = 50, 100 and 200. 
We compare the numerical light-curves with analytic results, and find that for a homogeneous 
external medium there is a simple scaling relationship among light-curves obtained for different 
parameters. The light-curves for power-law external densities are similar in shape to those in 
the homogeneous case. We discuss the implications of a comparison of the results with observed 
Gamma-Ray Burst time histories. 

Subject headings: gamma-rays: bursts - hydrodynamics - methods: numerical - relativity - shock 
waves 



Introduction 



Although the cosmological origin of Gamma-Ray Bursts (GRB) is still contested, there is good 
evidence to support it. The limits on the isotropy of their spatial distribution from the Burst and Transient 



Source Experiment (BATSE) put increasingly strong constraints on alternative hypotheses (Fishman & 



Meegan 1995] ; Meegan et al. 1996), the lack of evidence for repetition is becoming significant (Tegmark 
ct al. 1996), there is possible evidence for correlation with Abell rich clusters (Kolatt & Piran 1996) and 



cosmological time-dilation (Norris et al. 199?;; Norris et al. 1996), as well as spectral softening effects may 
have been found ( Mallozzi et al. 1995 , see however Mitrofanov et al. 1996 ). In this paper we proceed on the 
assumption that GRBs are cosmological in origin, and therefore result in a relativistically expanding fireball 



that produces a burst when shocks develop in the optically thin phase of the expansion (see Meszaros 1995 
for a review). Internal shocks occur at smaller radii, are marginally relativistic in the comoving frame, 
and can give rise to arbitrarily complicated light-curves whose duration is determined by the details of 



the energy release mechanisms (Rees & Meszaros 1994; Papathanassiou & Meszaros 1996). On the other 
hand, external shocks arising at larger radii from interaction with an external medium (EM) are extremely 
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relativistic in the forward blast wave that moves into EM (Meszaros, Rees & Papathanassiou 1994; Sari & 



Piran 1995). In this case, the burst duration is determined, for reasonable efficiencies, by the energetics 
of the source (E, rj) and the density of the external medium (n), via the Doppler time-delays of the light 
arriving from different parts of the light cone ( Rees 1966 ; Meszaros fc Rees 1993 ). 



Numerical simulations of relativistic jets with Lorentz factor ~ 22 have been published by Marti et 



ah (1995) and shock propagation over a limited region, generated from flows with Lorentz factors up to 
~ 2000 have been used as tests ( Marti fc Mullcr 1996 ; Romero et al. 1996 ), but numerical hydrodynamic 
simulations of shocks with very high (r < 50) Lorentz factors and their propagation over times that are 
several orders of magnitude longer than the system's crossing time, do not appear to have been published 
so far. This problem is of interest not only for GRB but also for relativistic pulsar winds interacting with a 
surrounding nebula, and for AGN jets, in which the intraday variability reported at radio wavelengths by 



Heidt fc Wagner (1996) , and at TeV energies by Gaidos et al. (1996) , may require bulk Lorentz factors in 



excess of 50-100. For GRB, such numerical calculations of the time variation of the dynamical variables are 
most urgently needed for computing accurate light-curves from external shocks. A discussion of BATSE 
3B pulse profiles has been presented by Fenimore, Madras & Nayakshin (1996), who model the burst as 
an energy release at constant power from a relativistically expanding surface. In this paper, we calculate 
numerically the energy release as a continuum of events from a relativistic fluid subject to shocks, assuming 
a radiating power that is proportional to the amount of internal energy stored in each elementary volume 
element. We discuss the dynamics of the forward and reverse shocks, including the phenomena of shock 
and rarefaction wave reflection, over the range of radii within which the ejecta randomize and radiate most 
of its initial kinetic energy. We compute numerical light-curves for initial bulk Lorentz factors r = 50 to 
200, in which the fireball ejecta interacts with a homogeneous or a power-law external medium, we compare 
these numerical light-curves with observations, and discuss the implications for GRB models. 



2. Initial Conditions in the Expanding Fireball 

The initial conditions of an impulsive fireball are: the total initial energy released E, the dimensionless 
entropy or initial random Lorentz factor rj = E/Mc 2 (or, equivalently, the amount of baryonic mass 
entrained in the fireball M = E/rjc 2 = 5.6 x 10 _6 £ ; 5i77^" M ) and the size ro of the region where the 
energy E is released. These parameters determine entirely the initial free expansion stage of the fireball 
(see Meszaros, Laguna & Rees 1993). Conversion of the energy E into radiation leads to a typical fluence 
at Earth of F = 10~ 6 E^D^ ergs/cm 2 , at a distance D — 10 28 Z?28 cm. In the beginning of the expansion, 
the fireball is optically thick with respect to pair formation and/or Compton scattering, out to radii as 
large as r ~ 10 13 cm in the source frame attached to the fixed origin point of the explosion (two other 
frames in use are the comoving frame attached to moving portions of the flow, and the detector or observer 
frame). As it expands, the fireball converts its internal energy into bulk kinetic energy. By the time it 
becomes optically thin, most of the initially available energy has been turned into kinetic energy, which can 
be reconverted into random energy and radiation via shocks. The above parameters and the characteristic 
particle density n of the external medium determine the deceleration radius around which this occurs, 
given by the condition E ~ (47r/3)r^ ec n77 2 rn p c 2 (see eq. ||). For computational economy, we do not start 
modeling the fireball-EM interaction from the very beginning of the fireball's expansion, but only after it 
has traveled for a few tenths of the deceleration radius. Before this stage is reached, the mass of the swept 
up EM is a small fraction of the mass within the deceleration radius; therefore any back-reaction caused by 
the interaction with the EM may be neglected and the fireball can be considered in free expansion. 
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During the free expansion phase most of the fireball's mass becomes concentrated in a thin shell near 
the outer edge of the expanding gas (Meszaros et al. 1993), and we will consider that practically all baryons 
are in this homogeneous shell, delimited by two discontinuities. The shell's Lorentz factor initially grows as 
r, and typically reaches a saturation value r ~ 77 at the saturation radius r s — nr Q — 10 8 772 7*0 g cm i after 
which the cooled fireball coasts at constant Lorentz factor. The lab frame ejecta shell width is constant 
A ~ ro until the broadening radius r b = rf r$ = 10 10 rj 2 ro,6 cm is reached, after which the shell's width 
increases linearly with radius: A ~ r/rj 2 . The optical depth for Thomson scattering r(r) = E n /Air r 2 c 2 77 
has the value r(r s ) = 3.5 x 10 10 E51 r^g r;^" 3 at r s (k = 0.40 cm 2 /g is the mass absorption coefficient 
for Thomson scattering) and is larger than 1 as long as log 772 + 0.67 logro,6 — 0.33 log £51 < 3.5 . 
The durations of the observed GRBs require that 77 ^ 10 3 , so that for ro ^ 5 x 10 9 cm and E = 10 51 
ergs, the Lorentz factor reaches saturation before optical thinness. Therefore the thinning radius is 
r t — (Ek/4itc 2 rj) 1 / 2 = 1.9 x 10 13 E^ n 2 1 ^ 2 cm. If r t is reached before r s , the shell's Lorentz factor at r t 
is less than 77, but subsequent Compton drag accelerate the shell up to a Lorentz factor comparable with 77 
(see Meszaros et al. 1993). From here on we will denote by V (~ 77) the maximum Lorentz factor the shell 
reaches at the end of the acceleration phase. Since the deceleration radius (see eq. S ) is much larger than 
both r s and r t , the shell travels at constant Lorentz factor T before the interaction with the EM cannot 
be neglected any longer and must be taken into account. For completeness, we note that rt is reached 
before r& if log n 2 + 0.4 logr ,6 — 0.2 log-E^i > 1.3 , and that the broadening radius is always larger than 
the saturation radius. 

The adiabatic index of the mixture protons + electrons + photons is determined by the component 
that gives the largest contribution to the total pressure. The ratio of the photon-to-electron (or proton) 
partial pressures is P v hoton I P electron = aVT 3 /3Nks, where N is the total number of electrons (baryons) 
and V is the comoving volume of the ejecta shell. Since the adiabatic index is initially 7 = 4/3 (the ejecta 
are radiation dominated, both in density and pressure), the product VT 3 is constant and the pressure ratio 
does not change. This means that the ejecta remain radiation dominated and the adiabatic index is 4/3 
as long as the electrons and photons are coupled, i.e. until the thinning radius r t is reached. After rt, 
the mixture electrons + protons is already cold (since r t > r s ) and the expansion is characterized by the 
adiabatic index 5/3. Taking into account the change in the adiabatic index after decoupling, the ejecta 
comoving density and pressure are given by: 

Po(r /r) 3 r<r s 

A)f/'~ 1 (?V r ) 2 r s <r<r b (1) 
PoV( r o/r) 3 r b <r 

Po(ro/ r ) 4 r<r s ( p {r n /r) 4 r < r s 

PoV~ 4/3 (ro/r) 8/3 r s <r <r t I Par]~ 4/3 (ro/r) 8/3 r s < r < r b 

jj ?7~ 4/3 (rt7^/r 5 ) 2 / 3 r t <r<r b P | p T] 4 / 3 {r a /r) 4 r b < r < r t ' 

Por) 2 (r 2 r^/r 15 ) 1 / 3 r b < r [ Po?? 4 ^ 3 (f , t?'o/ r5 ) r t < r 

if r t < r b and if rt > r b 

where po — 8.8 x 10 8 E51 r^g rj^ 1 g/cm 3 and po = 2.1 x 10 28 E§[ 4 r g 5 ^ 4 77^ 1 dyne/cm 2 are the baryonic 
density and pressure at r = ro. 

When the ejecta shell interacts with the EM, a forward shock (blast wave) sweeps up the EM, heating 
and collecting it in a dense sub-shell between the shock and the contact discontinuity (CD). A reverse shock 
propagates backwards into the fireball ejecta shell, compressing and heating the ejecta between the reverse 
shock and the CD. These shocks convert an increasing fraction of the ejecta bulk kinetic energy into internal 
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energy. The deceleration of the fireball is important when the energy stored in the shocked material is a 
substantial fraction of the fireball's initial kinetic energy (~ E). The random (or thermal) Lorentz factor of 
the shocked EM protons is comparable to the bulk Lorentz factor of the upstream ejecta (see eq. [3] and 
[10] in Blandford & McKee 1976). Therefore, for a homogeneous EM, the deceleration radius is given by 

r dec = {ZE/^jfnmpC 2 ) 1 ^ ~ 2.5 x 10 16 (S 5 i/n ) 1/3 % 2/3 cm , (3) 

where no is the EM number density in cm -3 . Since the shell is ultra-relativistic, the source-frame 
deceleration time-scale (or hydrodynamic time-scale) is given by t dec — fdec/c . Note that the 
deceleration radius is orders of magnitude larger than the saturation and thinning radii, so that 
at rdec the ejecta are very rarefied and optically thin. In the absence of the deceleration caused 
by the EM, the nominal lab frame shell width, comoving density and pressure at the deceleration 
radius, derived from equations (|lj), (0) and (||), would be: A d ec = 84 El( 3 tj^ 8 ^ 3 n^ 1 ^ 3 c s, 
Pdec" — 5-0 r\\ n c~ 2 ergs/cm 3 , and p s d h e f l — 1.5 x 10~ 15 E^^ 12 rj^ 2 r\\ nj 3 dynes/cm 2 if r t < r b 
or p s Jlf l — 1.9 x 10~ 14 E\[ 12 r^J nj 3 dynes/cm 2 if r t > rb- Throughout this paper, we will use the 
comoving rest mass energy density pc 2 , rather than the comoving mass density p, which has the advantage 
of showing whether the material is cold (p < pc 2 ) or hot (p > pc 2 ) . The width and the density and pressure 
inside the shell at the point from which we start modeling the interaction with the EM can be found from 
their nominal values at the deceleration radius, using the scaling equations p oc r -5 , p oc r -3 and A a r 
(from eqs. § and §). 

As viewed by the observer, the kinematics of the shell is similar to that of the relativistically expanding 
radio source described by Rees (1966). Due to the strong beaming effects, the observer receives radiation 
mainly from a spherical cap of opening angle ps r . This cap has an apparent speed towards the observer 
of ~ 2 r 2 c, so the effects of the deceleration of the ejecta shell are observed on a time-scale 

T burst = t dec I 2 T 2 = 42 El' 3 r~ 8/3 n 1/3 s. (4) 

Here we use lower case t to denote the lab frame (or source) time and upper case T for the detector time. 
Since the Lorentz factor of the shocked material will be less than the initial Lorentz factor T of the colliding 
shell, due to deceleration, the actual duration of the burst will be somewhat longer than predicted by 
equation ([|). This means that F should be in the range 10 2 — 10 3 in order to explain the observed burst 
durations. 



Numerical Results and Light Curves from Shocks with T 

External Medium 



100 in a Homogeneous 



Analytical approximations for the dynamics of highly relativistic shocks in the adiabatic limit and 



neglecting gradients of the physical quantities, were presented by Blandford & McKee (1976) and for GRB 
by Meszaros & Rees (1993), Meszaros et al. (1993), and by Sari & Piran (1995), Sari, Narayan & Piran 
(1996). The latter consider details of the dynamics for situations where the shocked material is hot and the 



reverse shock is either Newtonian or ultra-relativistic, as seen from the colliding shell's frame. Numerical 
results on relativistic shocks (mainly related to AGN jets or supernova remnants) have so far been restricted 
to cases with bulk Lorentz factors T & 5 — 10. Our goal is to model numerically the deceleration of the shell 
and the dynamics of the reverse shock for any value of its Lorentz factor, including the trans-relativistic case 
of considerable importance for bursts arising from internal shocks or from reverse shocks (in the external 
shock model). 
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For our numerical simulations we have developed a hybrid code based on standard Eulerian finite 
difference techniques in most of the computational domain, and a Glimm algorithm including an exact 
Riemann Solver in regions where discontinuities are present. The code was tested using standard shock 
tube problems and planar and spherical relativistic shock reflection problems. An analysis of the errors and 
convergence rates shows that the quality of our results is comparable to that in the Piecewise Parabolic 



Method (Marti & Miiller 1996; Romero et al. 199E). A description of this code and details of the tests 



performed to calibrate it are given in Wen, Panaitescu & Laguna (1996). The average mass and total energy 



relative errors during the runs done to obtain the numerical results described below, were less than 3%. 

In this section we outline the calculations for an ejecta shell colliding with a homogeneous EM, for a 
representative set of parameters T — 100, E — 10 51 ergs, ro = 10 8 cm, n = 1 cm~ 3 . Results for different 
parameters and the effect of an inhomogeneous EM are presented in §4. The evolution of the interaction 
shell-EM is followed from t = 0.4 td ec , at which point the physical parameters characterizing the shell are 
A = 33.5 cs = 1.01 x 10 12 cm, p = 5.74 x 10~ 12 dyne/cm 2 and p = 78.3 c~ 2 ergs/cm 3 = 8.70 x 10~ 20 g/cm 3 . 
Deceleration effects should be negligible up to this point since only 6.4 % of the EM mass within the 
deceleration radius has been swept up. We consider two simple cases: the adiabatic case (no energy leaves 
the system) and a simplified radiative collision case where it is assumed that a fraction 1/250 of the internal 
energy is radiated by the system every 8t — I0~ 3 tdec- We refer to this recipe for energy release as the 
"0.25 tdec energy release time-scale". 

Figure 1 compares the density, pressure and Lorentz factor 7 (we use the lower case letter for the flow 
Lorentz factor after the free expansion phase) inside the composite structure shocked shell-shocked EM, 
in the adiabatic and radiative cases. The abscissae give the (Eulerian) position relative to the CD that 
separates the ejecta and the shocked EM shells. Therefore, in all graphs, the CDs at different times between 
t = 0.5 t^c and t = 0.9 tdec, in steps of 0.1 tdec, are coincident. Each density profile has the same structure, 
showing from left to right: the unshocked fireball, the reverse shock propagating radially inward into the 
fireball, the shocked fireball material (condensed by a factor ^ 10), the CD between the inner (fireball) and 
outer (EM) shells, the shocked EM (« 100 times denser than the unshocked EM) and the blast wave (or 
forward shock) that propagates radially outward in the EM. The post-shock pressure and density satisfy the 



strong shock equations in Blandford & McKee (1976). One sees that pc 2 > p in the inner shell (therefore this 
shell is cold) and that p > pc 2 in the outer shell (the shocked EM is hot). If no energy leaves the structure, 
the internal energy of the shocked regions increases as more and more heated material accumulates between 
the two shocks, and accounts for the lost kinetic energy. Note that, in the adiabatic interaction, the shocked 
EM density is almost constant throughout the shell and slowly decreases in time, while in the non-adiabatic 
case the outer shell density is larger before the CD than behind the forward shock, and increases in time. In 
the adiabatic case, the thermal Lorentz factor of the shocked EM changes little with position and decreases 
from « 32 at t = 0.5 td ec to ~ 26 at t = 0.9 tdec- In the non-adiabatic case, the same thermal Lorentz factor 
is lower before the CD (« 26 at t — 0.5 tdec, ~ 10 at i = 0.9 tdec) than behind the forward shock (« 30 
at t — 0.5 tdec, ~ 20 at t = 0.9 tdec), since material shocked earlier had more time to radiate its internal 
energy. By the time the structure reaches 0.92 Tdec, the reverse shock has swept up all the ejecta shell gas 
(in both cases). The shock crossing time t^oss = 0-52 tdec — 4.4 x 10 5 s is within a factor 2 to the crossing 
time calculated by Sari & Piran (1995). The reverse shock velocity 0^ in the frame of the unshocked fluid 
(which moves in the lab frame at = Vr 2 — 1/r) is related to the lab frame reverse shock velocity 1 ^ by 
- 0[f ~ 02 r- 2 /(l - c ° h ), horn which c ° h = T 2 {0 - l s f)/[l + T 2 {0 - l s f)] (the speed of light is set to 
1). In the lab frame, - l 8 f = A/t s c ^ ~ 7.7 x 10~ 5 , therefore 0% = 0.43 (Lorentz factor 7™ ~ 1.10). 
The reverse shock is trans-relativistic (or quasi-Newtonian) in the ejecta shell rest frame. 
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After the reverse shock has crossed the fireball ejecta, a rarefaction wave propagates forward into the 
shocked inner shell. This wave travels in the comoving frame with the local sound speed, which in the 
adiabatic case is (3^ = c s = y/^p/h ~ 0.21, where h is the comoving enthalpy density. In the lab frame, the 
wave's velocity relative to the CD (which moves at 7cd ~ 66) is f3 l * h - j3 cd ~ + /%°) ^ 4.0 x 10" 5 . 

When the reverse shock crosses the inner shell, the width of this shell is Aj„ ~ 6.4 light-seconds, therefore 
the wave lab crossing time is ijfj."^ = A i „/(/3^f b — j3 c d) = 1.6 x 10 5 s ~ 0.19 tdec- In the non-adiabatic case, 
the wave's speed relative to the CD is j3 l ^ b — (3 c d ^ 3.7 x 10~ 5 and the inner shell's thickness is A in ~ 3.4 
light-seconds, therefore the wave lab crossing time is tj^oss = 9-0 x 10 4 s ~ 0.11 tdec- Figure 2 shows that 
after the rarefaction wave has crossed the inner shell, a second reverse shock forms and propagates in the 
now rarefied material behind the CD (this is easier to see in the Lorentz factor graphs, which show the 
shocked material being decelerated across the shock). Note that when energy is released from the system 
the shocked EM is denser. In the adiabatic interaction, the shocked EM behind the blast wave moves 
faster than that next to the CD, and a rarefaction wave develops in the outer shell. This is due to the 
fact that in the adiabatic interaction the shocked EM is hotter and a fraction of the internal energy is 
re-converted into kinetic, accelerating the forward shock more efficiently than in the non-adiabatic case. 
This effect can be seen in all graphs in figures 2 and 3, by comparing the forward shock position at the 
same times. In both situations (adiabatic or non-adiabatic), the thermal Lorentz factor in the outer shell 
decreases in time, due to the expansion of this shell (and the energy release, in the non-adiabatic case). 
Figure 3 shows the evolution of the structure until t = 2.0 tdec- At t = 1.92 tdec an the fireball's gas has 
been compressed behind the second reverse shock, and a new rarefaction wave propagates forward in the 
shell. At t = 2.0 tdec, when we stop the simulation, the whole structure extends over w 1000 light-seconds 
and consists of an extended rarefaction fan, a thin dense zone behind the CD and an extended outer shell. 
Note the much higher bulk Lorentz factor near the leading edge of the outer shell, in the adiabatic case. 

The light-curves generated during the fireball-EM collision were computed using a simplified energy 
release prescription which should give the main features expected from a more realistic radiation production 
treatment. We assumed that the lab-frame time-scale for radiating the internal energy is 0.25 tdec, in the 
sense defined above, and that radiation is emitted isotropically in the comoving frame. As shown by Rees 
(1996), the surfaces emitting photons that reach the observer simultaneously, are prolate ellipsoids with an 
axis ratio equal to the bulk Lorentz factor of expansion of the radiating medium. Due to strong beaming 
effects, only those regions moving at an angle less than T _1 relative to the axis origin-observer give a 
significant contribution to the light the detector receives. This corresponds to half of the surface of the 
above mentioned ellipsoids and, because the ellipsoids are very elongated, the light received by the observer 
at any given moment comes from a surface which can be extended over as much as one deceleration radius. 
Figure 4 (upper graph) shows the light-curve obtained by integrating the radiation emitted by a burst 1 
Gpc away from the observer, without taking into account cosmological effects. The inner curves represent 
the flux generated at times less than t = 0.5 tdec (the innermost curve), t = 0.6 tdec, t = 0.7 tdec and so on 
up to t = 2.0 tdec (the outermost curve), in steps of 0.1 tdec- This allows one to see when photons emitted 
at a certain time during the fireball-EM interaction arrive at the detector: light emitted before t = 1.2 tdec 
contributes to the GRB peak, while the energy released after t = 1.2 tdec forms the tail. In this case the 
burst lasts w 500 seconds and has a peak at T p ~ 80 seconds. Most of GRB's light (92%, see figure 4, 
middle graph) comes from the shocked EM (the main heat store). The log- log light-curve in figure 6 (upper 
graph) shows that the burst rise is a 1.2 power-law in the range 40 — 80 seconds, and that its tail (often 
referred to as "exponential decay") can be well approximated by a —1.2 power-law in the range 100 — 1000 
seconds. 

It is useful to compare these numerical results with the analytical calculations of Fcnimore et al. 
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(1996). They obtained impulsive GRB light-curves assuming that: (1) at any time, the radiating shell can 
be approximated as an infinitesimally thin surface, (2) the shell radiates at a constant comoving power Pq 
from to = 2 T 2 T until t max — 2 j 2 T max and is zero otherwise, and (3) the shell moves with a Lorentz factor 
that is a power-law in the detector time T : 7 = T (T/T )~^ . Fenimore et al. (1996) derived that the photon 
fluxes (photons/cm 2 s) should rise slower than T 2 and decay steeper than T~ 2 . The above assumptions are 
useful in order to develop an analytical model, but they imply substantial departures from the situation 
represented by our numerical results. In our case, the emitting structure is several hundreds light-seconds 
thick after t = 1.4 tdec, so that the light emitted by the gas moving toward the observer extends over the 
whole tail of the light-curve and, consequently, the structure cannot be approximated as infinitesimally 
thin. We do consider that radiation is emitted for a finite time (from to = 0.4 tdec to t max = 2.0 tdec), but 
at a power that is not constant: it rises, has a peak around t = 1.0 tdec and then falls (figure 5, middle 
graph). Moreover, different parts of the shocked structure move with different velocities. An averaged flow 
Lorentz factor, defined as the non-thermal energy-to-mass ratio, cannot be approximated as a power-law 
in T over the whole GRB history, but only locally. Nevertheless, one can try to ignore these complications 
and use the assumptions listed above to compute analytically the GRB bolometric flux (ergs/cm 2 s) and 
compare it with the numerical results. Similar to the derivations in Fenimore et al. (1996), we find that the 
energy flux at detector is: 



One difficulty arises from the fact that the shell is not very thin and T max , defined as the time when 
light emitted by the fluid moving directly toward the observer reaches the detector, is not well defined: 
radiation emitted at t = t max is spread in the light-curve over more than 500 seconds. If T is the time when 
light emitted from the forward shock arrives at the detector, then T max = 400 s, while if the CD is used 
to define T max , then T max = 910 s. Since the shell is very thin at to = 0.4 tdec, To has an unique value of 
To = 17 s. The power-law exponent £ for the average Lorentz factor is between 0.15 and 0.28 if 17 s < T < 
30 s, and between 0.53 and 0.61 if 100 s < T < 1000 s, depending on whether the blast wave or the CD is 
used to define T. Then, according to equation (||), the light-curve should rise as T a with 0.3 < a < 1.1 
for T £ 25 s, and should decay as T~ a with 1.2 < a < 1.8 for 300 s < T < 1000 s, and much steeper (as 
T~ 4 ) for T > 1000 s. Thus we obtained some analytical results that are close to the numerical ones (see the 
upper graph in figure 6). We emphasize the fact that the rise of the burst and the most important part of 
its tail are in the To < T < T max regime in equation (^|) , and that only the much steeper decay at T > 5000 
s in figure 6 corresponds to the T~ 4 fall for T > T max in equation 

The kinetic, internal and radiated energies, normalized to the total initial energy E = 10 51 ergs, and 
the radiated power in units of E/tdec = 1-19 x 10 45 ergs/s are shown in figure 5 (upper and middle graphs). 
For the 0.25 tdec energy release time-scale considered, 87 % of the initial kinetic energy was converted into 
internal energy until t = 2.0 tdec, from which more than 96% was released as radiation. As mentioned 
before, in this case the burst duration is determined by the how fast the kinetic energy is transformed 
into heat, i.e. the deceleration time-scale. Had we chosen a time-scale for energy release larger than the 
dynamical time-scale, the burst duration would have been determined by how fast the structure is able to 
radiate its internal energy, and the burst would last longer, having a lower fluence if a substantial fraction 
of the internal energy is lost adiabatically. 




T < Tc 



To < T < T, 
T < T 

x max _ - 1 







max 
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4. Discussion 

For given parameters E, r , r\ and n, the shape of the light-curves is basically independent of the 
details of the energy release mechanism, provided that the time-scale for radiating this energy is shorter 
than the deceleration time-scale (they should also be similar if the energy release time-scale exceeds the 
dynamic time-scale, but the two scale similarly with radius; otherwise the light-curves would be affected by 
an additional r-dependent efficiency factor tdyn/Uad)- For a homogeneous EM, the burst is smooth, has 
a sharp rise and a long power-law tail (rather than exponential). Figure 6 (top) is a log-log plot of the 
light-curves for E = 10 51 ergs, tq = 10 8 cm, n = 1 cm -3 and different initial Lorentz factor T (or parameter 
i]): 50, 100 and 200. It can be noticed that, as T is increased, the light-curve is shifted in the log F-log 
T plane, but its shape remains almost unchanged. Furthermore, the shift is the same for both pairs of 
successive parameters rj. The peak fluxes scale as 77 s / 3 and the times T p , T\ and T 2 at which the peak flux, 
and a fraction e~ x or e~ 2 of the peak flux are reached, scale as rj~ 8 / 3 (this is in fact what is expected from 
eq. g). A similar result is obtained if 77 and E are held constant and the EM density n is changed or if rj 
and n are fixed and E is varied (see figure 6, bottom graph). The peak fluxes scale as E 2 / 3 n 1 / 3 while T p , 
T\ and T2 go as E 1 ! 3 -nT 1 ! 3 . We conclude that, with a good accuracy, the light-curves can be parameterized 
in E, rj and n and reduced to a universal function /(r): 

F(E, r), n; T) — E 2 ^ 3 rf^ 3 n 1 ^ 3 /(r) , where t = E~ 1/3 ri 8/3 n 1/3 T. (6) 

This important feature of the external shock burst light-curves can be obtained analytically by substituting 
T (oc to /if cx t dec /i] 2 ) oc £1/3^-8/3^1/3 ( eq _ ||) in equat i on f or T < T < T max (which represents 
most of the light-curve rise and fall). Then the function /(r) is proportional to Po r]~ 2 nT 1 t 2 ~ 6 ^ , which 
implies that Po must be proportional to rj 2 n, in order to remove the multiplying factor dependence on burst 
parameters. This is indeed the case, since we assumed that the comoving radiated power is proportional 
to the comoving internal energy density of the shocked gas, which Blandford fc McKee (1976)| show to 



be proportional to rf n. Nevertheless, it must be emphasized that the assumptions used in the analytic 
derivations are rough, and the most reliable proof for the existence of such a simple scaling relation for 
F(T) is based on the numerical results. 

The effects of an inhomogeneous EM on the dynamics and energetics of the collision and on the 
light-curve can be estimated using a simple model of an expanding fireball interacting with an EM whose 
density varies as power-law (also of interest for the case where the event occurs inside a pre-ejected wind 
or nebula). Figure 4 (lower graph) shows the light-curves generated by the collision with an EM whose 
density is increasing as r 2 5 or decreasing as r" 1 , starting from r = OAr^ec (where n — lcm~ 3 ), compared 
with the homogeneous EM case. It can be seen that the light-curve is very insensitive to a monotonously 
changing EM density. The reason for this can be understood from figure 5 (middle and lower graphs), 
representing the shocked structure radiating power and average Lorentz factor as functions of the radial 
coordinate (the deceleration time-scale t dec and radius r,j ec are those previously defined for a homogeneous 
EM). When the external density increases with radius, the shocked gas on either side of the CD radiates 
more efficiently than in the homogeneous case, but the average Lorentz factor is lower and the emitted light 
is more stretched in detector time T than in the homogeneous EM case. The two effects act against each 
other, resulting in a light-curve that peaks at almost the same time and flux as for a homogeneous EM. A 
similar situation occurs when the external density decreases with r -1 (or with r -2 , as in a constant velocity 
wind): the two shocked regions radiate less efficiently, but the average Lorentz factor is larger, generating 
in the detector a flux with a very similar T-dependence. The impact of a fireball on an inhomogeneous, 
blobby medium ( Fcnimore ct al. 1996| ) would therefore lead to similarly shaped light-curves, with individual 



sub-pulses having the same generic shape as the light-curves in a homogeneous medium. We conclude that 
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bursts arising from external shocks have a fast rise and a long decay, and arc insensitive to variations in the 
EM density, as long as the energy release time-scale is shorter than the deceleration time-scale. 

I ~l 

Link & Epstein (1996) studied the temporal asymmetry of over 600 bursts from the BATSE 3B catalog, 

using the third moment of the burst time profile, and found that about two-thirds of them have positive 
time-asymmetry, in the sense that the flux rises faster than it falls. Ncmiroff et al. (1994)| used the brightest 



GRBs longer than 1 second detected by BATSE up to March 1993, and found that there are significantly 
more bursts that rise faster than the decay, on all energy bands. The average of the ratio fall-to-rise times 
over the full sample of bursts was found to increase with the time bin size, from 1.1 on the 64 ms time-scale 
to 1.4 on the 4096 ms time-scale. Mitrofanov et al. (1996) used 338 events of the Second BATSE catalog, 
that are longer than 1 second, aligned the normalized fluxes of all bursts at the peak time, and calculated 
the "averaged curve of emissivity" (ACE). Such an averaged light-curve represents the general time history 
of GRBs and is appropriate for a comparison with our numerical results (single hump light-curves) , since its 
construction averages out all individual peaks within bursts and allows the calculation of a time-asymmetry 
ratio that is not affected by the asymmetry of these peaks. Mitrofanov et al. (1996) computed the rise time 
Tfj and fall time Tp of the ACE by adding the flux weighted time bin durations in the ACE's rise and fall, 
respectively. They found that Tr = 2.33 s and Tp = 4.21s, which gives an asymmetry ratio Tp/Tr = 1.81. 
We note that the fall of the ACE (shown in figure 2 in their article) can be fit by a power-law of index -1.0, 
with time measured from T p — 1.87 s, where T p is the peak time. 

Using equation (Q) and the light-curve for T = 100 (figure 4, top graph), we find that the fall of the 
ACE published by Mitrofanov et al. (1996) can be best fit by the tail of the light-curve generated from an 
external shock with T w 400. The time-scales and asymmetry ratio of this burst are Tr = 0.90 s, Tp = 4.03 s 
and Tp/Tji = 4.50. Its fall is a -1.16 power-law in time measured from T p — 2.06 s. In conclusion, our model 
can reproduce well the tail of the ACE mentioned above, but the rise of the numerical light-curve is faster 
by a factor ~ 2.6, leading to a burst that is more asymmetric than the reported observations. Moreover, 
the external shock model with constant energy release time-scale always gives profiles with a decay longer 
than the rise, whereas a small fraction of the observed bursts shows the opposite trend. This implies either 
that additional factors are contributing to the light-curves from external bursts (e.g. perhaps a radiative 
precursor propagating ahead of the shock) , or a significant fraction of bursts arise from internal shocks where 



the light-curve is dictated by the details of an extended energy release that generated the fireball (Elees 



fc Meszaros 1994; Fenimore et al. 1996; Papathanassiou fc Meszaros 1996). The latter conclusion is also 



reinforced by the fact that internal shock bursts would be in principle able to show not only very different 
light-curve envelopes but also short time variability, whereas the external shock light-curves calculated here 
result in smooth profiles. This is due essentially to the smearing introduced by the simultaneous reception 
of radiation from regions moving at Lorentz factors 7 that differ by factors ~ 2. Nevertheless, our model 
should generate less asymmetric light-curves, and even ones with a rise front lasting longer than the tail, 
if specific energy release processes are taken into account and if these processes have time-scales that are 
longer than the dynamical time-scale in the beginning of the shell-EM collision (t < td ec ), and decrease as 
the interaction progresses, becoming shorter than tdec at later times (t > tdec)- 

In conclusion, we have developed a one-dimensional, finite difference -I- Riemann solver hydrodynamic 
code which is able to simulate the propagation of shocks with Lorentz factors of few hundreds and to evolve 
the shocked structure over a time that is ~ 10 4 F| times longer than its lab frame crossing time. We have 
used this code to simulate external shocks in cosmological Gamma-Ray Burst models arising from an initial 
impulsive energy release, and to analyze the dynamics of the forward and reverse shocks, as well as the 
shock and rarefaction wave reflections. These computations confirm the general conclusions from previous 
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simplified analytic descriptions of shock dynamics and light-curves, and provide detailed quantitative 
information which can be employed for comparisons with the observational data. We have also found a 
simple scaling law for bursts generated by external shocks with different physical parameters. 

This research has been supported in part through NASA NAG5-2362, NAG5-2857, NSF-PHY 93-09834 
and NSF-PHY 93-57219 
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Fig. 1. — Density, pressure and flow Lorentz factor for T = 100, E = 10 51 ergs, r = 10 8 cm and n = 1 cm~ 3 , 
at times indicated in the legend. The left column shows these profiles for the 0.25 td ec energy release time- 
scale, while the right column is for the adiabatic interaction. The structure is much thinner than its curvature 
radius and the position inside it is indicated relative to the contact discontinuity. Negative values correspond 
to the inner shell, positive values to the outer shell. Note that in the adiabatic interaction, the outer shell 
is less dense and more extended, and that the gradients in density, pressure and Lorentz factor are smaller. 
After t = 0.9 tdec the reverse shock crosses the inner shell, in both cases. 
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Fig. 2. — The same profiles in the energy release and adiabatic cases, after the reverse shock has crossed the 
fireball. The second reverse shock can be seen easier in the Lorentz factor 7 graph. 




Fig. 3. — Same as in figures 1 and 2, until t = 2.0 tdec, when most of the initial kinetic energy has been 
radiated (in the non-adiabatic case) and the simulation is ended. Shortly before t = 2.0 tdec, the second 
reverse shock crosses the rarefaction fan behind the contact discontinuity. 
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Fig. 4. — Upper graph: light-curve for T = 100 and 0.25 tdec energy release time-scale. The innermost 
curve is the light emitted by the shocked matter until t = 0.5 tdec, while the envelope is for the light 
emitted until t = 2.0 tdec- Between these two curves are shown the light-curves generated until t = 0.6 tdec, 
t = 0.7 tdec, t = 0.8 tdec in 0.1 tdec steps. The light-curve peaks at T v ~ 80 seconds. Middle graph: the 
shocked fireball and shocked EM contributions to the total light-curve. Note the large difference between the 
energies radiated by the two shocked media. Lower graph: comparison between the light-curves generated 
by a homogeneous EM (solid curve) and a power-law EM (dotted and dashed curves) . 
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Fig. 5. — Upper graph: Time evolution of kinetic (thick dotted line), internal (thin dotted line) and radiated 
energy (thin solid line), as the structure travels from 0.4 rdec to 2.0 rdec- All energies are normalized to the 
initial kinetic energy of the fireball. Middle graph: radiating power in units of E/tdec for homogeneous and 
power-law EM density, tdec and rdec as defined for a homogeneous EM are used, to allow comparison. The 
radiated power peaks earlier if the EM has an increasing density, as one would expect. Lower graph: average 
Lorentz factor for the same EM densities. E = 10 51 ergs, r = 10 s cm and r] = 100. 
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Fig. 6. — Upper graph: light-curves for E = 10 51 ergs, r$ = 10 s cm, n = 1 cm -3 and different initial Lorentz 
factors r (or parameters 77). Note that all light-curves have the same shape and shift equally when T is 
doubled. Two regions, were the received flux is a power-law in the detector time T, are identified. Lower 
graph: light-curves for 77 = 100, r = 10 8 cm and : E = 10 51 ergs, n = 1 cm -3 (thick solid curve); E = 10 51 
ergs, n — 0.01cm -3 (thin dotted curve); E = 10 51 ergs, n = 100cm -3 (thick dotted curve); E — 10 49 ergs, 
n = 1 cm -3 (thin dashed curve); E = 10 53 ergs, n = 1 cm -3 (thick dashed curve). 



